library(osmdata)
library(opentripplanner)
library(sf)
library(tidyverse)
library(ggthemes)
library(ggspatial)
library(viridis)
library(dplyr)

Introduction

I downloaded data points of all housing developments completed in 2020 that include at least 20 affordable units. In this exercise, I will investigate spatial relationships between these housing developments and stops along the T.

Tstops <- st_read("C:\\Users\\emmac\\Documents\\GitHub\\emmacolley-vis\\T_stops\\MBTA_Rapid_Transit.shp")
affordable_2020 <- st_read("C:\\Users\\emmac\\Documents\\GitHub\\emmacolley-vis\\affordable_2020\\massbuilds-shp-20201004-a44e92.shp") %>%
  filter(MUNICIPAL == "Boston") %>%
  st_transform(crs = "+proj=longlat +datum=WGS84") #project to wgs84


MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"
opq(bbox = 'Boston MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_xml(file = 'OTP/graphs/default/boston_streets.osm')
path_otp <- otp_dl_jar("OTP")
## The OTP will be saved to OTP/otp.jar
path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")

otp_build_graph(otp = path_otp, dir = path_data, memory = 1024)
## 2020-10-07 08:13:53 Basic checks completed, building graph, this may take a few minutes
## The graph will be saved to C:/Users/emmac/Documents/GitHub/emmacolley-vis/OTP
## 2020-10-07 08:14:07 Graph built
otp_setup(otp = path_otp, dir = path_data, memory = 1024)
## 2020-10-07 08:14:07 OTP is loading and may take a while to be useable
## Router http://localhost:8080/otp/routers/default exists
## 2020-10-07 08:15:08 OTP is ready to use Go to localhost:8080 in your browser to view the OTP
# Connect to opentripplanner
otpcon <- otp_connect()
## Router http://localhost:8080/otp/routers/default exists

1) Creating Isochrones

I plotted the new housing developments with affordable units and created two isochrones: a five minute walk and a five minute drive. I’ve also layered in the T stops (all train lines) to begin to see how close they are to the housing developments.

iso_5min_walk <- 
  otp_isochrone(otpcon = otpcon, fromPlace = affordable_2020,
                mode = "WALK", cutoffSec = 300) %>%
                st_transform(crs = MA_state_plane) %>%
  mutate(mode = "walk")
## Warning in otp_isochrone(otpcon = otpcon, fromPlace =
## affordable_2020, mode = "WALK", : Failed to get isochrone
## with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 300},"id":"fid--375ee30c_17502e3008f_-7fe8"}]}
iso_5min_drive <- 
  otp_isochrone(otpcon = otpcon, fromPlace = affordable_2020, 
                mode = "CAR", cutoffSec = 300) %>%
                st_transform(crs = MA_state_plane) %>%
  mutate(mode = "drive")

iso_all_modes <- rbind(iso_5min_drive, iso_5min_walk)
right_side <- st_bbox(iso_all_modes)$xmax
left_side  <- st_bbox(iso_all_modes)$xmin
top_side <- st_bbox(iso_all_modes)$ymax
bottom_side <- st_bbox(iso_all_modes)$ymin

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "cartolight", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = affordable_2020) +
  geom_sf(data = Tstops, color = "tomato1") +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_manual(values = c("skyblue", "royalblue"), name = "Area that is reachable\nwithin 5 min", labels = c("by car", "by walking")) +
  theme_map() +
  theme(legend.background = element_rect(fill = alpha("white", .5))) +
  labs(caption = "Red dots represent T stops.\nBasemap Copyright OpenStreetMap contributors.")
## Loading required namespace: raster

2) Plotting Walking and Driving distances

I’ve labeled each housing development with their neighborhood to begin to visual any trends between neighborhoods and commuting distances. Which neighborhoods have larger walk sheds? Which have large drive sheds? It’s hard to distinguish a pattern since I have a small sample set.

iso_areas <- iso_all_modes %>%
  mutate(area = st_area(iso_all_modes)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area) 

ggplot(iso_areas, 
       aes(x = as.numeric(walk), y = as.numeric(drive))) +
  geom_point(color = "lightslategrey", size = 3) +
  geom_text(aes(label = c("Roxbury", "South Boston", "South Boston", "Dorchester", "Mission Hill", "Brighton", "Jamaica Plain","Jamaica Plain","South End","South End")), hjust=-.2, size = 3) +
  scale_x_continuous(name = "Area within a five-minute walking distance\nof housing development (square km)",
            breaks = breaks <- seq(10000, 200000, by = 20000),
            labels = breaks / 1000000) +
  scale_y_continuous(name = 
            "Area within a five-minute driving distance\nof housing developmen (square km)",
            breaks = breaks <- seq(0, 700000, by = 100000),
            labels = breaks / 1000000) +
  theme_minimal()
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).

3) Assigning a Transit Score to Each Isochrone

I will use now investigate how many T stops are within each isochrone to see which housing developments are more accessible to the T by walking or driving. I am including both isochrones in this study.

st_as_sf(iso_all_modes)
iso_transformed <- st_transform(iso_all_modes, crs = "+proj=longlat +datum=WGS84")
iso_transformed <- iso_transformed %>%
  mutate(transit_score = lengths(st_covers(geometry, Tstops)))
ggplot(iso_transformed) +
  annotation_map_tile(zoomin = 1, type = "cartolight", progress = "none") +
  geom_sf(aes(fill=transit_score), alpha=.5) +
  scale_fill_viridis_c(option = "C", name = "Walk and drive sheds\nwith T stops", breaks = breaks <- seq(0, 4, by = 1), labels = paste(prettyNum(breaks), "T stops")) +
  theme_map() 
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Ok, that figure is a bit helpful but it looks like most of the development’s walk and drive sheds only have 1 T stop within them. I will take a step back and compare isochrones with and without T stops.

4) Which Isochrones overlap with a T stop?

6 of the 10 isochrones have a T stop within them. That means it would take 5 minutes to walk or 5 minutes to drive to the T from 6 of these affordable housing developments. They are shaded in blue in the figure below.

iso_transformed <- iso_transformed %>%
  mutate(num_Tstops = lengths(st_covers(iso_transformed, Tstops))) %>%
  mutate(has_Tstops = num_Tstops > 0)
## although coordinates are longitude/latitude, st_covers assumes that they are planar
n_Tstops_iso_transformed <- sum(iso_transformed$has_Tstops)

n_Tstops_iso_transformed
## [1] 6
ggplot(iso_transformed) +
  annotation_map_tile(zoomin = 1, type = "cartolight", progress = "none") +
  geom_sf(data = iso_transformed,
          aes(fill = has_Tstops)) +
  scale_fill_manual(values = c("gray85", "darkblue"),
          name = "Boston Neighborhoods\nby presence of a community center", 
          labels = c("Without a T stop",
                     "With a T stop")) +
  theme_map()
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

otp_stop()